The Molecular Hubbard Hamiltonian: Field Regimes and Molecular Species 
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The molecular Hubbard Hamiltonian (MHH) naturally arises for ultracold polar alkali dimer 
molecules in optical lattices. We show that, unlike ultracold atoms, different molecules display 
different many-body phases due to intrinsic variances in molecular structure. We also demonstrate 
a wide variety of experimental controls on molecules via external fields, including applied static 
electric and magnetic fields, an AC microwave field, and the polarization and strength of optical 
lattice beams. We provide explicit numerical calculations of the parameters of the MHH, including 
tunneling and direct and exchange dipole-dipole interaction energies, for the molecules ^Li^^^Cs, 
2^Na^°K, ^'^Rb^^^Cs, ^^K^'^Rb, and ^Li^^Na in weak and strong applied electric fields. As case 
studies of many-body physics, we use infinite-size matrix product state methods to explore the 
quantum phase transitions from the superfluid phase to half-filled and third-filled crystalline phases 
in one dimension. 



Molecules cooled to ultralow temperatures provide ac- 
cess to a wealth of exciting new physics, such as the study 
of chemical reactions at ultracold temperatures, preci- 
sion tests of fundamental symmetries, novel means of 
quantum information processing, and many-body physics 
with many internal degrees of freedom [1 . The vast tun- 
ability proffered by molecules has led to many propos- 
als to manipulate molecules with external fields in order 
to simulate outstanding unsolved quantum Hamiltonians 
motivated by solid state physics [IHi] . Taking a comple- 
mentary approach, we address the immediately accessible 
many-body physics observable in upcoming experiments 
using ultracold molecules in optical lattices, in which ex- 
ternal fields will not be fine-tuned. In this Letter, we 
use the molecular Hubbard Hamiltonian (MHH) to show 
that the first experiments with ultracold molecules in op- 
tical lattices will display a variety of many-body phases 
which are strongly dependent on molecular species and 
amenable to control with external fields. 

For ultracold alkali atoms, interactions are well- 
modeled by a contact pseudopotential proportional to 
the s-wave scattering length, and the 5-wave scattering 
length is readily tunable by using a Feshbach resonance. 
Hence, while the specific positions of Feshbach resonances 
may differ among different alkali species, essentially the 
same tunability of interactions is provided for all species. 
We demonstrate that the same is not true for ultracold 
molecules. Due to intrinsic variations in the structure 
of alkali dimer molecules, different molecules access dif- 
ferent regimes of the MHH, and so display a variety of 
many-body phases. In this work, we take as examples 
the bosonic molecule ^^Rb^^^Cs [5] and the fermionic 
molecules ^Li^^^Cs [6^, ^SNa^^^j^ [7, ^^K^'^Rb ^, and 
^Li^^Na [9 , all of which are being explored in current 
experiments. After identifying the MHH for a particu- 
lar near-term experimental setup we discuss some of its 
experimental controls which are not available to ultra- 
cold atoms. We then provide explicit numerical values 
for the parameters appearing in a quasi one-dimensional 



reduction of the MHH such as tunneling and dipolar in- 
teraction terms. To provide case studies of many-body 
physics, we numerically calculate quantum phase tran- 
sitions from superfluid to crystalline phases at rational 
densities in one dimension using variational matrix prod- 
uct state techniques. By computing phase diagrams in 
terms of experimental parameters such as the strength of 
an applied electric field, we show that crystalline phases 
will not be seen for the molecules with the weakest dipole 
moments. 

The molecular Hubbard Hamiltonian. The molecular 
Hubbard Hamiltonian (MHH) in the simplest experimen- 
tal configuration is 

Here a labels the internal state of a molecule, and the op- 
erator CLi^cr {a\^) represents either a hard-core bosonic or 
fermionic annihilation (creation) operator of a molecule 
in state a on lattice site i. In order, the terms in Eq. ([T]) 
are tunneling of molecules between sites i and j with 
state-dependent tunneling energy tj-i^cr^ local transitions 
between molecular internal states a ^ driven by an 
AC microwave field Eac which has circular polarization 
in the lab frame and effective Rabi frequency ftaa'^ en- 
ergetic detuning A^r of a molecule in state a from a 
reference ground state, direct dipole-dipole interactions 
Uj-i^cra' between a molecule in state a on site i and a 
molecule in state a' on site j, and exchange dipole-dipole 
interactions Ej_i^(j^^^^>^^'_^ which exchange a quantum of 
rotation between molecules at sites i and j through the 
dipole- allowed transitions ai and cf2. The 

parameters r^, r^/, and rE represent cutoffs of the range 
of tunneling and direct and exchange dipole-dipole inter- 
actions, respectively, and set a consistent energetic order 
of approximation for the terms in Eq. ([T]) for a given set 
of field parameters. We note that these cutoffs generally 
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depend on the internal state, rt = n,cr etc., and on direc- 
tionality, e.g. the range of tunneling will be shorter along 
directions in which the lattice potential is deeper. 
The Hubbard parameters in Eq. ([T]) are defined as 

tj-i,- = - /bz dqe-'^-C^'-'^^^E^JvBz , (2) 
A<T = -io,<T - w<5yv,i , (3) 
^ (U^^^ + 5f /) ± (U^X + , (4) 

1^^^/ = Excdi^aa' , = / c^r^J^ (r)4^0a' (r) , (6) 

where vbz is the volume of the first Brillouin zone (BZ), 
Ri is the coordinate of lattice site i, and the positive (neg- 
ative) sign in Eqs. Q and ([5| refers to bosons (fermions). 
The Wannier functions ^^;^c^(^) are the quasimomentum 
Fourier transforms of the lowest band Bloch functions 

^spV^aq(r) = £^aqV^aq(r) , (7) 

which are the simultaneous eigenfunctions of the lattice 
translation operators and the single-particle Hamiltonian 
^^sp- The dipole and short-range interaction integrals 
appearing in Q and ([5| are 

(8) 

^^^A ^ J (r)P (a) ff}, (r') , (9) 

where /^*^/(r) = w^^{T)wi'a' {'^) is a Wannier function 
product, dq is the projection of the dipole operator along 
space-fixed spherical basis element e^, Cq (R) is an un- 
normalized spherical harmonic, and P{a) is a regularized 
5- (p-) wave pseudopotential [10] for bosons (fermions) 
which is assumed to have the universal dependence on 
the van der Waals length a discussed in Ref. j^l \l2\ • We 
evaluate the integrals Eq. ([8| using the ab initio method 
of Ref. [13]. As shown in that work, the dipole-dipole 
interaction parameters Uj-i^aa' and ^j_^^cricr2cr^cri have a 
^ l/\j — i\P power-law tail with p ^ 3 and a possibly 
significant exponentially decaying contribution when the 
lattice strengths are anisotropic. 

The single particle Hamiltonian consists of kinetic, lat- 
tice, and internal components, and the internal contribu- 
tion has been discussed in detail in past work [M] . We use 
the few- molecule parameters in Table [l] which have been 
calculated from density functional theory [T^ or avail- 
able experimental data The internal states of an 
alkali dimer molecule in the presence of a DC electric 
field Edc = £'00^2 group into manifolds Mn^v). 
Here, N denotes that the manifold adiabatically corre- 
lates with having quanta of rotation as £^dc 0, Mn 
is the projection of rotation along e^, and rj indexes the 
nuclear hyperfine degeneracy. Two such manifolds with 
I Mat I 7^ I^atI are well separated in energy provided that 



TABLE I: The few- molecule parameters can differ by sev- 
eral orders of magnitude for experimentally relevant 
molecules 1121 [TB, 16 . In order, we give the permanent dipole 
moment d, the nuclear spins Ii and I2, the rotational con- 
stant Bn, the nuclear quadrupole coupling constants {eQq)i 
and {eQq)2, the nuclear ^-factors gi and ^2, the hyperfine spin 
couplings ci, C2, C3, and C4, the magneto-association field 
and the van der Waals length a. 
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the scaled DC field strength /3dc = dE^c/BN > 0.2. 
In addition to Edc, we consider an applied DC mag- 
netic field B = Be^. The values of B we use, collected 
in Table |l] correspond to the positions of the Feshbach 
resonances used for magneto-association [5 -9 ; these are 
the magnetic fields naturally and easily used in exper- 
iments. The Zeeman splittings and intrinsic couplings 
between hyperfine sublevels produce intra-manifold split- 
tings which are typically a few tens to hundreds of kHz. 

In order to induce rotational transitions between states 
in manifolds |0,0,7^) and |1,±1,7^'), we require a single 
AC microwave field with frequency uj and circular po- 
larization with respect to the DC field axis; such fields 
are already used in experiments ^ . The detuning A^r 
denotes how far from a single-molecule resonance state 
a is, and the effective Rabi frequency Q^a' expresses the 
strength of transitions from state a to state a' . In strong 
B fields where the intra-manifold splitting is typically 
larger than e.g. the dipole-dipole interaction energy, the 
AC field determines which states are dynamically acces- 
sible. Hence, the number of states involved in the MHH 
dynamics can be modified by the static field configuration 
and the frequency and power of the AC field, resulting 
in a tunable quantum complex system [14 . The num- 
ber of dynamically accessible internal states can be quite 
sizable, typically a few tens, due to a large hyperfine de- 
generacy, see Ii and I2 in Table [ij 

Finally, we consider that the molecules are subject 
to a simple cubic optical lattice potential consisting of 
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FIG. 1: (Color online) Examples of experimental controls not 
available to atoms, (a) An applied electric field ^dc gives 
rise to single-particle states with tunable expected dipole mo- 
ments. Vertical dashed lines indicate the position of Edc = 
5kV/cm for the various molecules, (b) The difference in tun- 
neling energy between internal states may be tuned by a vari- 
ety of controls which have no counterparts for atoms such as 
the DC field strength ^dc , the angle between the lattice po- 
larization e and Edc, and the frequency u of the microwave 
field. 

three pairs of counter-propagating laser beams, Matt(r) = 
^r]e{x,y,z} cos^(7r7//a)eJ; • a • e^. Here, a is the polar- 
izability tensor of the molecules, E'^ is the intensity of 
the laser field along direction 77, is the polarization of 
the field along 77, and a = 545nm is the lattice spacing. 
In the numerical examples we choose the polarizations 

= = = corresponding to current exper- 

iments [17]. The ratio of the components of the polar- 
izability due to coupling to S and 11 states at the given 
wavelength, /a^^ is taken to be 0.223 [18]. Due to the 
fact that the molecular polarizability tensor is rank two, 
non-e^ polarized light couples together the states in the 
manifold |1,±1,77) [19 and the tunneling amplitude tp^cr 
depends on the internal state and the polarization, see 
Table [n] and Fig. [ijb). Also as a result of the coupling in 
the manifold |1, ±1, 77), the single-particle eigenstates are 
in general a quasimomentum-dependent superposition of 
the states which diagonalize the internal Hamiltonian in 
free space, and so the dipole moments generally cannot 
be taken outside the integration in Eq. ([8|. 

Tunability of MHH parameters for two internal states. 
The given experimental setup provides several controls 
which have no counterpart in atoms. To provide explicit 
examples of this tunability, we now focus on a special case 
of the MHH in which only two internal states are popu- 
lated, a e {1,2}, and the lattice is strongly confining in 
the y and z directions. The condition of only two internal 
states is well-realized for current experiments with large 
B fields and weak microwave power [51 [16]. The state 
|cr = 1) is taken to be the absolute ground state, and 
the state |cr = 2) is the lowest energy eigenstate which 
maximizes the di transition dipole moment with \a = 1) 

m- 

Figure [TJa) shows the dependence of the dipole mo- 
ments dq^aa'i scc Eq. ([6|, ou the scaled DC field strength 
/3dc- The resonant moments {do^cra}^ which are the 
expected dipole moments of the states {|cr)}, give rise 
to the direct dipole-dipole interactions Uj-i^aa'^ and 



TABLE II: Order of magnitude differences in Hubbard pa- 
rameters give rise to a variety of phases. Here we collect the 
many-body parameters of the MHH Eq. ([T|) for Vy = Vz = 
55Er, Vx = SEr, the few- molecule parameters given in Ta- 
ble^ and the ^dc and species given in the table. The value 
^DC = 5kV/cm is a typical strong field in experiments. We 
do not include Q, as it is widely tunable by adjusting the 
power of the applied AC microwave field. 
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the transition moments {di^crcr'l^which give the strength 
of a dipole- allowed transition a ^ a', yield the ex- 
change terms ^j_^^cricr2cr^cri • Thc transition moments are 
strongly affected by hyperfine structure and lattice po- 
larization effects as shown by the discrepancy between 
<^i,i2 and di^2^ latter being computed in the absence 
of a lattice. In Fig. [ijb) we display the tunneling ener- 
gies of two internal states as a function of the quasi- ID 
lattice depth Vx = {cr = l\Ele^^ ■ a ■ ex\(J = 1) measured 
in units of the recoil energy Er = t:'^ / 2mo? . The dis- 
parity between the tunneling in different internal states 
can be altered by changing the strength of the DC field, 
tilting the lattice polarization, or changing the frequency 
of the AC field, as indicated by the various curves. 

Many-body case studies. To demonstrate the emergent 
many-body features of Eq. ([T]) in near-term experimen- 
tal setups, we present many-body case studies employ- 
ing the infinite-size variational ground state search using 
matrix product states (IMPS) for ID systems [21]. We 
choose the energetic cutoff defining etc. in Eq. ([T]) to 
be IHz, and the field configurations are chosen such that 
terms not included in Eq. ([T]) such as assisted tunnel- 
ing and tunneling outside of the quasi- ID domain fall 
below this threshold. These cutoffs account for the fact 
that many-body phases induced by long-range interac- 
tions may be fragile against other terms in the Hamilto- 
nian, non- conservative effects such as heating from the 
optical lattice, and lifetime constraints of the ultracold 
gas. 

We begin by considering the phase diagram of Eq. 
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in a quasi- ID geometry and absent an AC field. Without 
the AC field only the absolute ground state a = 1 is 
populated, and so the MHH reduces to 

H=- Eij.r, h-M^j + I Eij;ru Uj-ifliflj . (10) 

At half filling, p = 1/2, a phase transition occurs at 
strong coupling (small to a crystalline phase (CP) 

with long-range order in the density correlation function 
A/'(r) = (AnoAn^), An^ = fii — p; an exponential de- 
cay of the single-particle density matrix A{r) = (alar)', 
and a non- vanishing single-particle gap corresponding to 
fixed density for a range of chemical potential /i, see 
Fig. [2] To within numerical resolution, we find that the 
CP melts directly into a superfiuid with algebraic decay 
of both JV{r) and A{r) and no single-particle gap, and 
there is no region of supersolid phase displaying simul- 
taneous long-range order in A/'(r) and algebraic decay 
of A{r). Due to their small dipole moment, ^^K^^Rb 
and ^Li^^Na do not display the CP phase transition for 
fields £^DC ^ 5kV/cm, and so are not included in Fig. |2j 
The absence of the CP transition for these molecules is 
an example of how different molecules display different 
many-body physics. 

When the interaction Ux is a convex function of x, 
it has been conjectured that regions of rational density 
other than p = 1/2 can occupy finite regions in the phase 
diagram [22 . To study the accessibility of such phases, 
we investigate the CP at filling p = 1/3 for ^Li^^^Cs in 
Fig. [2|f). Here we fix £^dc and vary the quasi-lD lattice 
height Vx in order to change the relative strengths of in- 
teraction and tunneling. While it is possible to achieve 
this phase with the other molecules, the equilibration 
timescales will be longer than the lifetime of typical ex- 
periments due to the smaller dipole interaction energy 
which sets the timescale, see Table [nl Hence, ^Li^^^Cs 
is the best choice for investigating phases due to strong 
direct interactions. Similarly, other CPs with densities 
p < 1/3 could be seen in principle, but the requirement 
of very weak tunneling with respect to interactions is 
suggestive that such configurations may be difficult to 
equilibrate on reasonable timescales as well as obtaining 
a sufficiently cold system. We note that the p = 1/3 CP 
is unique to molecules with long-range dipole-dipole in- 
teractions; there is no analog of this phase in atoms with 
short-range interactions. 

The simplest situation utilizing the internal structure 
of the molecule is when there is exactly one molecule per 
lattice site and only two internal states are accessible. 
Due to the hard-core constraint the translational degrees 
of freedom are frozen out and only the internal degrees of 
freedom are dynamical. When the two internal states of 
the molecule are mapped to an effective two-component 
spin via the standard transformation = {hu — n2i)/2, 
S^' = a|^a2i, S~ = {S^y , the MHH for either bosons or 
fermions maps, up to constant terms, onto a long-range 




FIG. 2: (Color online) Phase diagram of the MHH in the 
absence of an AC field. (a) Density phase diagram for 
LiCs, showing the p — 1/2 crystalline phase (CP) with non- 
vanishing single-particle gap. The dashed line is a guide to 
the eye for the region of fixed p = 1/2. (b,c) analogs of (a) for 
NaK and RbCs, respectively, (d) Density correlation func- 
tion J\f{r) with period-two oscillation removed for the points 
marked in (a), (e) Single-particle density matrix A{r) for the 
points marked in (a), (f) Density phase diagram for LiCs near 
the p = 1/3 CP transition at fixed ^dc = 3.5kV/cm showing 
small density fluctuations induced by our use of the grand 
canonical ensemble [21]. 

XXZ model in transverse and longitudinal fields 

H =i Eij,r. x^-Sts- - 2n sf + sf (ii) 

Here, Ij^i = Uj-i^u - 2Uj-i^i2 + ^j-i,22, Xj-i = 
£^i-i,i2i2, and hz = (Ai - A2) + Ej-i;r^ - 
Uj -1^22)- For the given experimental setup, both the 
Ising-type coupling and the XX- type coupling Xj_i 
are positive, corresponding to antiferromagnetic interac- 
tions. This model has been studied in ID recently with- 
out the Vt term [23]. As is increased relative to the 
other terms in the Hamiltonian, the system transitions 
to a trivial paramagnetic phase with nonzero magnetiza- 
tion along X, corresponding to strong on-site correlation 
of the internal state. We stress that the ratios of the 
various spin couplings depend in an essential way on the 
molecular species, see Table [Hj 

We have presented the molecular Hubbard Hamil- 
tonian, the natural Hamiltonian governing the many- 
body physics of ultracold molecules in optical lattices. 
Molecules are sensitive to experimental controls which 
are not relevant for atoms such as an applied DC elec- 
tric field, the polarization of optical lattice light, and the 
frequency and power of an AC microwave field driving 
rotational transitions. The response of a molecule to ex- 
ternal fields is also much more dependent on molecular 
species than is the case for atoms. Our results establish 
that the molecules employed in current experiments will 
display a range of species-dependent many-body features. 
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Supplementary material to "The Molecular Hubbard Hamiltonian: Field Regimes and 

Molecular Species" 

M. L. Wall, Erman Bekaroglu, and Lincoln D. Carr 

INFINITE-SIZE VARIATIONAL GROUND STATE SEARCH USING MATRIX PRODUCT STATES 

Infinite-size variational ground state search using matrix product states (IMPS) is a numerical method which finds 
a representation for the ground state of an infinite one-dimensional lattice model |1J. The state is assumed to be 
invariant under shifts by a user- specified number of lattice sites and hence can be represented in terms of a 
periodically repeating unit cell with q lattice sites. The unit cell is decomposed as a matrix product state (MPS) 
with finite entanglement cutoff x nieasured by the Schmidt rank [2 , and this unit cell is iteratively optimized via the 
solution of an eigenvalue problem. The optimization reaches an end when the density matrix p obtained by tracing 
over all sites to the right of a given point in the infinite lattice changes by less than a tolerance e between successive 
iterations. We measure that change via the orthogonality fidelity 



^ortho = Tr Y V PnPn-lV Pn , (12) 

where the subscript on p denotes the optimization cycle index. In practice, we choose the orthogonality fidelity 
tolerance to be e = eiocab where eiocai is unit-cell averaged 2-norm distance between the optimal unit cell and its 
representation as an MPS with fixed entanglement cutoff x- That is, deviations of the density matrix p between 
successive iterations are only from the errors due to a finite cutoff of entanglement. Once the unit cell has converged, 
a standard orthonormalization procedure is used to define a normalized state which conforms to a known canonical 
MPS form on the infinite lattice [3]. From this state, any desired observable may be computed. 

After optimizations have been performed with a range of entanglement cutoffs x, the effects of finite entanglement 
cutoff can be discerned by a scaling analysis in x- It is known rigorously that MPSs with a finite cutoff x cannot 
display power-law decay of correlation functions for arbitrary distances. Hence, for a given x, there exists a distance 

beyond which a correlation function with a power-law decay begins to decay exponentially, and is an increasing 
function of x- In addition, the scaling for some universal quantities, such as the correlation length, have a known 
scaling relation with x which depends only on the data of the conformal field theory describing criticality in the 
infrared limit |4]. An example of the scaling of correlation functions with x is given in Fig.lsjwhich shows the scaling 
of the single-particle density matrix A{r) in the superfiuid and crystalline phases of Eg. (|1Q[) near half filling. All 
results given in the main text have been checked for convergence in x and the unit cell length. 

1 
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FIG. 3: Scaling of single-particle density matrix with entanglement cutoff in iMPS. The single-particle density matrix A{r) 
computed by iMPS for x =12 (red), 16 (green), 20 (blue), and 40 (purple) in the crystalline and superfiuid phases. In the 
crystalline phase, A{r) decays exponentially and so the results converge for finite x- In the superfiuid phase A(r) has an 
algebraic decay, and the range over which this decay is well- approximated increases with x- 

Long-range interactions are accommodated by representing the Hamiltonian as a matrix product operator 
(MPO) [5 . When the interaction cutoff ru is less than 6, we use a finite-range interaction on ru sites, but for 
larger cutoffs we use an infinite-range Hamiltonian obtained by fitting the function Up to a series of exponentials [6] . 
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We do this out of numerical convenience; the infinite-range decomposition requires fewer numerical resources for the 
same level of accuracy when r^/ > 6. 

The density plot near the p = 1/3 crystalline phase, Fig. 2(f) of the main text, appears noisy. This apparent 
noise is an artifact of having a finite resolution of points in the plot and using the grand canonical ensemble. We 
demonstrate this point in Fig. [4j Panel (a) displays the density near the phase transition, and panel (b) shows the 
same plot except that the region to the left of the solid vertical line has been computed with double the resolution in 
chemical potential. The density fluctuation spots become smaller for higher resolution, and are concentrated on the 
border between two regions of fixed density. Hence, as the resolution increases, the effect of these fluctuations will 
become less prominent. Density fluctuations on the border between superfluid regions of fixed density are also seen 
for finite systems in the grand canonical ensemble, but do not appear in the canonical ensemble 0. Finite-size effects 
introduce a no n- vanishing single-particle gap, and this inherently limits the resolution of the chemical potential for 
the grand canonical ensemble in a finite- volume system. 




fi/Ui ii/Ui 

FIG. 4: Number fluctuations in the grand canonical ensemble, (a) Density plot of the region near p = 1/3. (b) Same as in (a), 
but the region to the left of the solid vertical line is computed with twice the resolution in chemical potential. 



[1] 1. P. McCulloch, arXiv:0804.2509; M. L. Wall, Ph.D. thesis, Colorado School of Mines, 2012. 

[2] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003). 

[3] U. Schollwock, Ann. Phys. (NY) 326, 96 (2011). 

[4] F. Pollmann et a/., Phys. Rev. Lett. 102, 255701 (2009). 

[5] L P. McCulloch, J. Stat. Mech., P10014 (2007). 

[6] M. L. Wall and L. D. Carr, New J. Phys., in press. 

[7] L. D. Carr et al Phys. Rev. A 81, 013613 (2010). 



